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ABSTRACT 

We build a stellar-dynamical model of the Milky Way barred bulge and disk, 
using a newly implemented adaptive particle method. The underlying mass 
model has been previously shown to match the Galactic near-infrared surface 
brightness as well as gas-kinematic observations. Here we show that the new 
stellar-dynamical model also matches the observed stellar kinematics in several 
bulge fields, and that its distribution of microlensing event timescales reproduces 
the observed timescale distribution of the MACHO experiment with a reasonable 
stellar mass function. The model is therefore an excellent basis for further studies 
of the Milky Way. We also predict the observational consequences of this mass 
function for parallax shifted events. 

Subject headings: Galaxy: kinematics and dynamics — Galaxy: bulge — Galaxy: 
disk — Galaxy: structure 



1. Introduction 

It is now known, from several independent lines of evidence, that the Milky Way Galaxy 
(MWG) is barred (e.g., Gerhard (2001)). However, a comprehensive model consistent with 
the main observables - luminosity distribution, stellar-kinematics, gas-kinematics, and mi- 
crolensing - has so far been still missing. Recently, Bissantz & Gerhard (2002) obtained 
a luminosity density model for the MWG from the dust-corrected L-band COi?i?/DIRBE 
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map of Spergel, et al. (1995), through a non-parametric constrained maximum hkehhood 
deprojection. This model (hereafter: COBE-p model) is consistent also with the observed 
magnitude distributions of clump giant stars towards several bulge fields, and with the mi- 
crolensing optical depth towards the bulge derived from these stars (Popowski, et al. 2003; 
Afonso, et al. 2003); see also Binney, et al. (2000) and Bissantz & Gerhard (2002). Further- 
more Bissantz, et al. (2003) found that the hydrodynamical gas flow in the potential of the 
COBE-p model matches the observed gas dynamics of the inner MWG well. 

The structure of the inner MWG can also be constrained by observations of stellar 
kinematics along fixed lines of sight (Sharpies, ct al. (1990) [Sh90]; Spacnhaucr, et al. (1992) 
[Sp92]; Minniti, et al. (1992) [Mi92]) and by the microlensing event timescale distribution 
(ETD) (Alcock, et al. 2000). The ETD has been studied largely with models which assume 
some distribution of disk and bulge kinematics {e.g. Han & Gould (1996); Peale (1998); Mera, 
et al. (1998)). An exception was Zhao, et al. (1996), who used the dynamical bar model of 
Zhao (1996) augmented by an analytic disk model, but failed to match the long duration 
(t > 100 days) tail of the ETD. In the present Letter we show that a full stellar-dynamical 
model based on the COBE-p model is consistent with these independent data as well. 

Dynamical models of the MWG have been generated using the Schwarzschild (1979) 
method, in which the distribution function (DP) of a galaxy is built from numerically inte- 
grated stellar orbits. Pollowing earlier work by Zhao (1996), Hafner, et al. (2000) constructed 
a 22168-orbit dynamical model of the MWG. Dynamical models of the MWG have also been 
obtained by A^-body methods (Pux 1997). Sycr & Trcmainc (1996) [ST96] introduced a novel 
method for generating self-consistent dynamical models. The Syer-Trcmaine (hereafter ST) 
method is allied to the Schwarzschild method, but, rather than superposing time-averaged 
observables from an orbit library, the ST method constructs a model by actively varying the 
weights of individual particles (orbits) as a function of time. This permits arbitrary geom- 
etry and a larger number of orbits to be used in the model building. Our dynamical model 
for the COBE-p density in the MWG is constructed with the ST method, demonstrating its 
usefulness for real galaxy modeling. This Letter compares the model's bulge kinematics and 
microlensing ETD with their observed counterparts. 



2. The Syer-Tremaine method 

The idea of the ST algorithm is to assign individual weights to particles of a simulation, 
which are then changed to reduce the deviation between the model and observations. An 
observable Yj associated to a stellar system characterized by a distribution function /(z), 
z = (x, v), can be written as Yj = J Kj{z)f{z)(fz, where Kj{z) is a known kernel. If this 
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stellar system is simulated with N particles having weights Wi and phase space coordinates 
Zi, then we can write the observables of the simulation as yj{t) — J^^^ Wi{t)Kj{zi{t)). ST96 
define the "force of change" on the weights as 



The small and positive parameter e governs how rapidly the weights are pushed such that 
the simulation observables yj{t) converge towards the observables Yj. The constants Zj act 
as normalizations. The full ST method also includes prescriptions for temporal smoothing 
and a maximum entropy term, to reduce fluctuations. We have implemented the ST method 
with the MWG disk-plane surface density as the observable (Debattista et al. in preparation). 
We set e — 0.25, a — 0.524, // = 0.001, where a and // are the parameters of the temporal 
smoothing and the entropy terms, respectively, in the notation of ST96. 



Since the MWG contains a bar, our initial model also had to be barred. The simplest way 
to achieve this was to evolve an A^-body model of an initially axisymmetric, bar-unstable, 
disk galaxy. The A^-body simulation which produced the barred model consisted of five 
disk and bulge components inside a frozen halo. The frozen halo was represented by a 
cored logarithmic potential. The initially axisymmetric disk was modeled by a truncated 
exponential disk. Disk kinematics were set up using the cpicyclic approximation to give 
Toomrc Q = 1.3. The disk and bulge were represented by 4 x 10® equal-mass particles, with 
a mass ratio : Mb = 4:1. Further details of the setup methods and model units can 
be found in Debattista (2003) [D03]. We use the halo, disk, and bulge parameters given in 
Table 2 of D03, which give a flat rotation curve out to large radii. 

The simulation was run on a 3-D cylindrical polar grid code (described in Sellwood & 
Valluri (1997)), with technical parameters exactly as in D03. The initially axisymmetric 
system was unstable and formed a rapidly rotating bar at t ~ 50. By t — 160, the bar 
instability had run its course and further secular evolution of the bar was mild. The resulting 
system did not match the COBE-p model of the MWG and needed to be evolved further 
with the ST code. First, however, we eightfold symmetrized the COBE-p model in order to 
reduce the amplitude of spirals, which we did not try to reproduce. We evolved the A^-body 
model from t = 160 under the ST prescription with the fixed potential of the COBE-p model 
plus dark matter halo. We kept the bar pattern speed at its value in the A^-body model, 
which scales to 56 km s~^ kpc~^, consistent with the MWG (Dehnen (2000); Debattista, et 
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al. (2002); Bissantz, et al. (2003)). At t = 240 {i.e. ~ 4 bar rotations), we shut off the ST 
algorithm and evolved the system to t — 280 to assure that the particles are phase mixed. 



To compare our dynamical COBE model (the COBE-Djn model) with observations, 
we adopted the same viewing parameters as were used to determine the COBE-p model: 
Rq — 8kpc, Zq — 14pc and </7bar = 20° (Bissantz & Gerhard 2002). We scale the velocities in 

the COBE-Djn model to the MWG by matching to the local circular velocity. We assumed 
that the local standard of rest has only a circular motion, with vlsr = 200 km s~^, and we 
adopted the values of the solar peculiar motion from Dehnen & Binney (1998). 

The densities of the COBE-Dyn and COBE-p models match very well, with azimuthally 
averaged errors smaller than 5% out to Rq. The largest errors (< 15%) occur in small isolated 
regions on the bar major-axis. In the (unconstrained) vertical direction, the disk is somewhat 
thicker than the MWG at Rq, but this leads to a change in optical depth r towards Baade's 
window of < 15%. In the bulge region, on the other hand, the scale-height of the COBE- 
Dyn model matches that of the MWG very well. We compared the model's kinematics to 
observations towards Baade's window (Sh90, Sp92) and in the field at (/, b) = (8°, 7°) (Mi92), 
using the selection functions determined by Hafner, et al. (2000). Table 1 shows our results. 
The overall fit of our model to the observed kinematics is rather good. 



We now show that the COBE-Dyn model is also consistent with the microlensing ETD. 
Alcock, et al. (2000) presented an ETD, corrected for their experimental detection efficiency, 
based on 99 events in 8 fields. Popowski (2002) argued that one of these fields seems biased 
towards long-duration events, introducing some uncertainty in the observed ETD. Here we 
use the full-sample Alcock et al. ETD in order that our results may be compared with 
previous ones. We computed the ETD with the self-consistent kinematics of the COBE-Dyn 
model. A microlensing event is characterized by the source distance, Ds, lens distance, 
Dl, the proper motion, v±, of the lens with respect to the line-of-sight between observer and 
source, and lens mass, M^. The probability P{i) for observing an event duration i = 2Qe/v± 
is given by 



3. Results: Density and Bulge Kinematics 



4. The Microlensing ETD 
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HML)v±fMS{i - 2QE/vA_)dv^ dDL dDs dML. (2) 

Here p{d) is the density of tiie MWG at distance d from the observer along the Hne-of-sight to 
the observed field, $(Mi) is the mass function (MF) of the lens population, Qe{Ds, -Dl, Ml) 
is the Einstein angle, and f{v±) is the distribution of v±. We solved the multiple integral by 
Monte Carlo random drawings of the parameters {Ds, Dl, v±, Ml) as follows: (1) To obtain 
the source distance (0 < Ds < D™'^^ = 12kpc), we used the COBE-p model, since this is 
less noisy than the particle realization. The probability of is cx 

to account for a magnitude cut-off (Kiraga & Paczyhski 1994). (2) The lens distance 
(0 < -Dj^ < Dg) was selected from p{Dl) j^^ p{D L)dD l, where p is a normalized probability 
density distributed as in the COBE-p model. (3) For the relative velocity v^, we used the 
particle distribution of the COBE-Dyn model, randomly selecting a particle at ~ Ds and 
another at ~ Dl- The proper motions of these particles then determined v±. (4) The lens 
mass Ml/Mq was selected from a Kroupa (1995) MF, $(Ml) = (3{Ml/Mq)-^, with 

( (2.35, 0.1038) M< < Ml/Mq < 0.35 
(7,/?) = J (0.6,0.6529) 0.35 < Ml/Mq < 0.6 (3) 
[ (2.35, 0.2674) 0.6 < Ml/Mq < M> 

We explored varying M^ and Mf". We obtained the ETD, shown in Fig. 1, by simulating 
10^ events, and weighting each by the remaining factors in Eqn. 2. We tested our Monte 
Carlo integrations by reproducing one of Peale (1998) 's model ETDs. 

We started with {M^, M^) = (0.075, 10), for which we obtained a Kolmogorov-Smirnov 
distance between data and model of Dks = 0.213. (We excluded the bin ait < 3.1 days from 
the MACHO data in all such comparisons, because it appears to be too heavily affected by its 
large detection-efficiency correction.) To improve on this fit, we first explored the effects of 
uncertainties in the COBE-p model. The most important of these is c/^bar- Setting (/9bar = 30°, 
we found only a minor change to the ETD, in agreement with Peale (1998). Making the 
bar stronger, or the disk velocity dispersion outside the bar smaller did not alter the ETD 
substantially. Therefore we next explored variations in the MF. Like Peale (1998), we found 
that modest changes can improve the fit substantially. Our best fit, with Dks = 0.068 was 
obtained with M^ — 0.04 and Mf" = 10. However, a more conservative limit is M£" = 4, 
which gives Dks — 0.081. (If the suggestion of Popowski (2002) is correct, which would shift 
the ETD peak to smaller i, then a smaller Mf" would be required anyway.) 

We now explore the causes of long duration (LD) events in the COBE-Djn model; using 
(M^, M^) = (0.04, 4) as our standard model for this analysis. We start by noting, from Fig. 
1, that the vast majority of sources are located in the bulge (6 < Ds < 10 kpc). This is 
also true for the lenses responsible for short duration events, but disk lenses become more 
important at longer durations; indeed, for i > 25 days, one third of the lenses are at D^, < 4 
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kpc. In Fig. 2 we separate the ETD into the near and distant lens sub-samples and show the 
heliocentric angular velocities and cumulative distributions of Mi for both. Note first that 
lenses with Mi > 0.5Mq contribute significantly to LD events in both the near and distant 
sub-samples. Lens mass, however, is not the full explanation of the LD events, as has been 
noted by previous studies, and the relative motions of lens and source in the heliocentric 
frame must also be considered. The kinematics of the LD sources are substantially those of 
a rotating triaxial bulge/bar which points almost towards the observer: thus their apparent 
tangential motions are largely due to the solar motion, giving fttan,s ~ 205/8 km s~^ kpc~^. 
Distant lens LD events are then possible because the lenses share very similar kinematics 
with the sources (note that massive lenses become necessary only in the last quartile, i > 60 
days). For the nearby lens sample, LD events have a rather large spread in fltan,L (due both 
to their proximity and the velocity dispersion of the COBE-Djn model), which together with 
larger M^'s is able to produce LD events. We conclude, therefore, that there is no single 
cause for the long-duration events. 

Standard 3-parameter fits to microlcnsing events are symmetric about the time of peak 
amplification, resulting in a degeneracy between Mi, v±, Di and Ds- One degree of de- 
generacy is removed by also measuring the light-curve shift due to the parallax from earth's 
orbit, which gives a relation between v±_ and Di/Ds- These shifts are present in all events, 
but most go undetected because of infrequent sampling and photometric errors. Buchalter 
& Kamionkowski (1997) estimate a 1% detection efficiency of parallax-shifted events for the 
MACifO-type setup, and much higher for second generation experiments. The fight curves 
of such events require 5 parameters, including k = R^{D2^ — Dg^)/QE, where i?® = 1 AU. 
In Fig. 3, wc present our predictions for the probability distribution in the (k, t)-plane, as- 
suming 100% detection efficiency. These distributions are twin-peaked, with the lower peak 
increasingly separated from the global peak as decreases, as it must since k oc 0^^ while 
t (xQe- The location of the second peak may therefore provide an observational constraint 
on the MF at low mass. 



5. Conclusions 

We have presented a dynamical model of the MWG constructed using the Syer-Tremaine 
method, constrained only by the MWG density map of Bissantz & Gerhard (2002). Although 
no kinematic constraints were used, the model (i) matches observed bulge kinematics in sev- 
eral fields and is (ii) able to reproduce the observed microlensing event timescale distribution. 
For the best fitting MF, the model (iii) predicts a twin-peaked probabifity distribution in 
the (k, i)-plane, which may be observationally tested with new generations of microlensing 



- 7- 



experiments. (iv) The underlying mass model has been previously shown to match the Galac- 
tic near-infrared surface brightness as well as gas-kinematic observations. It is therefore an 
excellent basis for further studies of the Milky Way 

This work was supported by the Schweizerischer Nationalfonds through grant 20-64856.01. 
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(8°, 7°) 


Mi92 


45± 10 


46 




(8°, 7°) 


Mi92 


85 ±7 


100 



Table 1: Comparison of kinematic quantities computed from the COBE-Dyn model with 
observations. In the first column, the superscript ^ indicates that the value given is helio- 
centric, and ^ that it is Galactoccntric. All quantities are in km except for a^^ and a^^, 
which are in milliarcseconds year~^. 
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Figure Captions 

Fig. 1 The ETD of the COBE-Djn model compared to the detection-efhciency-corrected 
observations of the MACHO group (histograms in all panels). Top: cumulative distribution 
function for the standard model, {M^,M^) — (0.04,4), (sohd line), the best model with 
{M<,M>) = (0.04,10) (dotted), and a model with {M<,M>) = (0.075,10) (dashed). We 
obtain D^s = 0.081, 0.068, 0.213, respectively, for the 3 models. Middle panel: differential 
distributions of these models (same line styles). (In these two panels, all model distributions 
have been smoothed with a kernel density estimator of bandwidth 0.1.) Bottom: ETD of 
the {M^,M^) = (0.04,4) model and its decomposition into events with: 6 < Ds < 10 kpc 
(dotted line), < Dl < 4 kpc (dashed) and Dl > 4 kpc (dot-dashed). 

Fig. 2 Top two panels: The long duration events {i > 42 days) in the (M^, M£") = (0.04,4) 
model. Both are maps (on the same relative scale) in the plane of heliocentric tangential 
angular velocities, fltan,L and fltan,s- Left: near lenses {Dl < i kpc); right: distant ones 
{Dl > 4 kpc). The diagonal and horizontal dashed lines indicate fltan,L = ^tan,s and 
^tan,s — 205/8 = 25.6 km s~^ kpc~^, respectively. Bottom: distributions of Ml for distant 
(solid lines) and nearby (dashed) lenses. The different hues result from splitting into quartiles 
by contribution to the full ETD the distribution of events sorted on t. Event durations 
increase as the mean mass increases. The heavy curve shows the underlying mass function. 

Fig. 3 Predicted probabihty distribution of parallax shifted events in the (k, f)-plane, for the 
standard model with (M^,M£") = (0.04,4). We use a smoothing kernel with (5K,(^iogt) — 
(0.01,0.1). The stars mark the locations of secondary peaks when — 0.075, 0.04, 0.02 
and 0.01 in order of increasing k. 
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Fig. 2.- 
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